Compressibility of solid helium 
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The compressibility of solid helium ("^He and *He) in the hep and fee phases has been studied by 
path-integral Monte Carlo. Simulations were carried out in both canonical (NVT) and isothermal- 
isobaric {NPT) ensembles at temperatures between 10 and 300 K, showing consistent results in 
both ensembles. For pressures between 4 and 10 GPa, the bulk modulus B is found to decrease by 
about 10%, when temperature increases from the low-temperature limit to the melting temperature. 
The isotopic effect on the bulk modulus of helium crystals has been quantified in a wide range of 
parameters. At 25 K and pressures on the order of 1 GPa, the relative difference between '^He and 
*He amounts to about 2%. The thermal expansion has been also quantified from results obtained 
in both NPT and NVT simulations. 
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I. INTRODUCTION 

In the last decades there has been a continuous 
progress in the study of different types of substances un- 
der extreme conditions of pressure and temperature, thus 
enlarging appreciably the experimentally accessible re- 
gion of phase diagrams pL]. In particular, the influence of 
controlled hydrostatic pressure on structural, thermody- 
namic, and electronic properties of various kinds of solids 
has been deeply studied. Pressures on the order of tens 
of GPa can now be routinely applied to real materials 
iiii. 

Solid helium, in spite of having been studied for many 
years, has a broad interest in condensed matter physics 
because of its peculiar character of "quantum solid" . In 
particular, its zero-point vibrational energy and associ- 
ated anharmonic effects are markedly larger than in most 
known solids [6]. In addition, its electronic simplicity al- 
lows one to carry out detailed studies, that would be 
enormously difficult for other materials 0, H, Q • The in- 
terest on the behavior of solids under high pressures has 
been also focused on solid helium. Thus, diamond-anvil- 
cell and shock-wave experiments have allowed to study 
the equation of state (EOS) of sohd ^He up to pressures 
on the order of 50 GPa |^, [lOJ, [ill. In last years, the 
effect of pressure on heavier rare-gas solids has also been 
of interest for both experimentalists [J, II4I and theorists 

MMMM- 

The Feynm an p ath-integral formulation of statistical 
mechanics [13, [lH is well suited to study thermodynamic 
properties of solids at temperatures lower than their De- 
bye temperature &£>, where the quantum character of 
the atomic nuclei becomes important. In particular, the 
combination of path integrals with computer simulation 
methods, such as Monte Carlo or molecular dynamics, 
has revealed as a powerful technique to carry out quanti- 
tative and nonperturbative studies of many-body quan- 
tum systems at finite temperatures. This has allowed to 
study several properties of solids further than the usual 
harmonic or quasiharmonic approximations |7[. 



The path-integral Monte Carlo (PIMC) method has 
been used to study several properties of solid helium 
[3, [ig . \2(X [21I [22 . [23I [24I a s well as heavier rare-gas 
solids [23, [2a [23, [2a, [23, [20]. For helium, in particu- 
lar, this method has predicted kinetic-energy values [1^| 
and Debye- Waller factors [31| i n g ood agreement with 
data derived from experiments [32, [33| . PIMC simula- 
tions have been also employed to study the isotopic shift 
in the helium melting curve |20l . l2l| . The EOS of solid 
helium at T = has been studied by diffusion Monte 
Carlo in a wide density range \§, ^54| , as well as at finite 
temperatures by using PIMC simulations with several in- 
teratomic potentials J22j |. 

In last years, there has been a debate on the existence 
of supersolidity in ^He at temperatures lower than 1 K 
[33, [SOj [33, [331 • This debate is still open, but is out of the 
scope of this paper, since we consider here solid helium at 
temperatures higher than 10 K, where quantum exchange 
effects between atomic nuclei are not relevant. 

In this paper, we study the compressibility of solid 
■^He and ^He by PIMC simulations. We employ the 
isothermal-isobaric (NPT) ensemble, which allows us to 
consider properties of these solids along well-defined iso- 
bars. For comparison, we present also results of PIMC 
simulations in the canonical (NVT) ensemble. By com- 
paring results for ^He and "'He, we analyse isotopic effects 
on the compressibility as a function of pressure. 

The paper is organized as follows. In Sec. II, the com- 
putational method is described. In Sec. Ill we present 
and discuss the results, that are given in several subsec- 
tions dealing with thermodynamic consistency, pressure 
and temperature dependence of the bulk modulus, iso- 
topic effects, and thermal expansion. Finally, in Sec. IV 
we present the conclusions. 



II. METHOD 

Equilibrium properties of solid ^He and ''He in the 
face-centred cubic (fee) and hexagonal close-packed (hep) 



phases have been calculated by PIMC simulations. The 
PIMC method relies on an isomorphism between the 
quantum system under consideration and a fictitious clas- 
sical one, obtained by replacing each quantum particle by 
a cyclic chain of Q classical particles {Q: Trotter num- 
ber) , connected by harmonic springs with a temperature- 
dependent force constant. This isomorphism appears as 
a consequence of discretizing the density matrix along 
cyclic paths, which is usual in the path-integral formu- 
lation of statistical mechanics [13, [lM|- Details on this 
computational method are given elsewhere 0, [S^, So, [4l| . 

Helium atoms were considered as quantum particles in- 
teracting through an effective interatomic potential, com- 
posed of a two-body and a three-body part. For the two- 
body interaction, we employed the potential developed 
by Aziz et al. ^ (the so-called HFD-B3-FCI1 poten- 
tial). For the three-body part we took a Bruch-McGee- 
type potential [43, [4J|, with the parameters given by 
Loubeyre [44| , but with parameter A in the attractive ex- 
change interaction rescaled by a factor 2/3, as suggested 
in Ref. [21| . This interatomic potential was found earlier 
to describe well the vibrational energy and cquation-of- 
state of solid helium in a broad range of pressures and 
temperatures (23l |. 

Our simulations were based on the so-called "primi- 
tive" form of PIMC [H, El]. We considered explicitly 
two- and three-body terms in the simulations. The ac- 
tual consideration of three-body terms did not allow us 
to use effective forms for the density matrix, developed 
to simplify the calculation when only two-body terms are 
explicitly considered [2l[ . Quantum exchange effects be- 
tween atomic nuclei were not taken into account, because 
they are negligible for solid helium at the temperatures 
and pressures studied here. (This is expected to be valid 
assuming the absence of vacancies and for temperatures 
higher than the exchange frequency ~ 10~^ K [3|.) To 
calculate the energy we have used the "crude" estimator, 
as defined in Refs. (39l.l45j. 

We have employed both the canonical [NVT) ensem- 
ble and the isothermal-isobaric {NPT) ensemble. Our 
simulations were performed on supercells of the fee and 
hep unit cells, including 500 and 432 helium atoms re- 
spectively. To check the convergence of our results with 
system size, we carried out some simulations for other su- 
percell sizes, and found that finite-size effects for N > 400 
atoms are negligible for the quantities studied here. In 
particular, changes of the bulk modulus with the size of 
the simulation cell were found to be smaller than the 
statistical error bar of the values derived from our simu- 
lations. 

Sampling of the configuration space was carried out 
by the Metropolis method at temperatures between 12 K 
and the melting temperature at each considered pressure. 
For given temperature and pressure, a typical run con- 
sisted of 10^ Monte Carlo steps for system equilibration, 
followed by 4 X 10^ steps for the calculation of ensem- 
ble average properties. Each step included attempts to 
move every replica of every atom in the simulation cell. 



In the NPT simulations, it also included an attempt to 
change the volume. To keep roughly constant the ac- 
curacy of the computed quantities at different tempera- 
tures, we took a Trotter number Q proportional to the 
inverse temperature, so that QT — 3000 K. This means 
that for solid helium at T = 20 K we had Q — 150, and 
a PIMC simulation for N = 500 atoms effectively in- 
cludes 75000 "classical" particles. More technical details 
are given in Refs. [23, [2^, [4g| . 

In the NVT ensemble the pressure and isothermal 
bulk modulus B = —V{dP/dV)T can be obtained from 
thermal averages of various quantities obtained in PIMC 
simulations, as shown in Appendix A. The derivation 
is straightforward from the partition function Z^vt for 
N quantum particles in the canonical ensemble, but is 
rather tedious in the case of the bulk modulus, due to the 
number of terms appearing in the volume derivatives. In 
particular, B can be obtained from Eqs. (|A11[) and (|A12[) 
in Appendix A. In the NPT ensemble, the isothermal 
bulk modulus is related to the mean-square fluctuations 
of the volume V of the simulation cell by the expression: 
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as shown in Appendix B. 



III. RESULTS AND DISCUSSION 
A. Consistency checks 

We have calculated the bulk modulus in the NPT and 
NVT ensembles. In general, we prefer the isothermal- 
isobaric ensemble, so that we can study solids along 
well-defined isotherms. NVT simulations are, however, 
frequently used in PIMC simulations [3, [l^, [3l|, and a 
comparison of results obtained in both ensembles seems 
necessary as a consistency check of the method, and in 
particular for the case of solid helium. 

We discuss first results obtained in the NPT ensem- 
ble. At constant pressure, the relative fiuctuations in the 
volume, ay /V , can be found from Eq. ([I]): 
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For hep ^He at 25 K we found in PIMC simulations 
av/V = 2.9 X IQ-^ for P = 2 GPa, and a lower value 
of 1.4 X lO^'' for 20 GPa (for a simulation cell including 
A^ ~ 432 atoms). This means that the product BV in- 
creases as the pressure is raised [see Eq. ([2], in spite of 
the reduction of V , indicating that the growth of bulk 
modulus with the pressure dominates in the product BV 
(see below). 

For a cubic crystal, the fluctuations in the lattice pa- 
rameter a can be derived from Eq. ([T]) to give: 
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FIG. 1: Relative fluctuation of the lattice parameter, aa/a, 
of fee "^He vs L~^'^ , for simulation cells of lateral size La. 
Shown are results of PIMC simulations for two pressures (2 
and 7 GPa) at T = 100 K, and for L = 3, 4, 5, and 7. Error 
bars are less than the symbol size. 
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FIG. 2: Isothermal bulk modulus of solid *He as a function 
of pressure at three temperatures: 25 K (squares), 150 K 
(circles), and 300 K (triangles). Open and filled symbols cor- 
respond to hep and fee helium, respectively. Error bars are 
smaller than the symbol size. Dashed lines are fits to the data 
points using the expression: B = Bo + BgP + ^B'^P^. 



where L^ is the number of cubic unit cells in a simulation 
cell with side length La. From Eq. ([3]) one can see that 
the relative fluctuation in the lattice parameter, aa/a., 
scales as L~^'^. This normalized fluctuation for fee '*He 
is displayed in Fig. 1 as a function of L~^^^. Shown are 
results of PIMC simulations in the NPT ensemble at T 
= 100 K and two pressures: 2 and 7 GPa. The linear 
dependence shown in this figure agrees with the depen- 
dence of a a on the simulation-cell size given by Eq. ^. 
The different slopes of these lines are basically due to the 
change of compressibility with pressure. 

To check the thermodynamic consistency of our simu- 
lations in both NPT and NVT ensembles we have car- 
ried out several tests. The first obvious test consists 
in taking the crystal volume derived from isothermal- 
isobaric simulations, and using it as an input in NVT 
simulations at the same temperature. The latter should 
give the same pressure [using Eq. (jA7|) ] as that used ear- 
lier in the NPT simulations. As an example, for fee 
^He at 25 K and 1 GPa, we find first a lattice parame- 
ter a = 3.6955 A, that introduced as an input in NVT 
simulations yields a pressure of 1.0002(5) GPa, in good 
agreement with the input in the previous NPT simula- 
tions. Going on with the same set of parameters, we can 
also check the consistency of the bulk modulus B derived 
from both types of simulations. At constant pressure {P 
= 1 GPa), we find Bnpt = 4.50(2) GPa from the vol- 
ume fluctuations [see Eq. (jB5p in Appendix B] vs Bnvt 
= 4.47(2) GPa derived at constant volume by using ex- 
pressions (IA11|) and (|A12p . 



B. Pressure and temperature dependence 

We now consider the pressure dependence of the bulk 
modulus, as derived from our PIMC simulations. This is 
shown in Fig. 2 for '*He at three temperatures. In this 
figure, T decreases from top to bottom: T = 25, 150, and 
300 K; open and filled symbols correspond to hep and fee 
helium, respectively. At each considered temperature, we 
present data for the pressure region where the solid was 
stable (or metastable) along the PIMC simulations, a re- 
gion that becomes broader as the temperature is lowered. 
In addition to the typical rise of B for increasing pres- 
sure, we find that for a given pressure, B decreases as 
the temperature is raised. In the pressure range shown 
in Fig. 2 one also observes a departure of linearity in the 
dependence of B on pressure. By fitting these results to 
the expression B — Bo + B'qP + ^BqP^, we can find the 
pressure derivatives at P = (Sg and Bq). Thus, for T 
= 25 K we obtain Bq = 0.47 GPa, B'^ = 4.01, and S^' = 
-0.068 GPa-i. 

For the other temperatures shown in Fig. 2 (T = 150 
and 300 K), an extrapolation of the results to P = 
yields negative values of Bq, indicating that the solid is 
mechanically unstable at these temperatures for low pres- 
sures. At 25 K and zero pressure, even though the liquid 
is known to be the stable phase, the solid can still be 
metastable, as it has not yet reached the limit of mechan- 
ical stability, or the spinodal line . Apart from this, the 
curves B{P) at the considered temperatures are rather 
parallel one to the other in the common stability region. 

It is interesting to compare the bulk modulus obtained 
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FIG. 3: Isothermal bulk modulus of solid ''He as a function of 
temperature at three hydrostatic pressures: 4 GPa (squares), 
7 GPa (circles), and 10 GPa (triangles). Open and filled sym- 
bols correspond to hep and fee helium, respectively. Error 
bars of the simulation results are on the order of the symbol 
size. Dashed lines are guides to the eye. 
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FIG. 4: Isothermal bulk modulus of solid ^He versus density. 
Open squares represent results of path-integral Monte Carlo 
simulations in the NPT ensemble. Error bars are less than 
the symbol size. A solid line indicates the bulk modulus de- 
rived by Zha et al. ^31 from Brillouin scattering experiments. 
The dashed line is a guide to the eye. 



here for solid helium with those of heavier rare-gas solids 
at the same conditions. For example, at T = 20 K and 
P = 1 GPa, we find for sohd ^He: B = 4.53 GPa, vs 7.20 
and 9.43 for Ne and Ar, respectively [SO]. This is in line 
with the known result that the bulk modulus increases 
with atomic mass for given values of temperature and 
pressure [SOJ . 

Rare-gas solids have been studied in metastable con- 
ditions, even at negative pressures 47[. They have been 
found to be metastable in PIMC simulations close to the 
spinodal point, defined as the point at which the com- 
pressibility K — 1/B diverges (or the bulk modulus van- 
ishes). Thus, for solid Ne and Ar at 5 K, the spinodal 
point was found at P — —91 and —245 MPa. For solid 
helium, we could not approach the spinodal point at any 
temperature, due to the large quantum fluctuations that 
make the solid to be unstable along the simulations. This 
happened even at zero pressure and relatively low tem- 
perature, in spite of the fact that the expected bulk mod- 
ulus at these conditions is still far from vanishing. For 
example, at 25 K we find for solid helium Bq = 0.47 GPa, 
as derived from extrapolation of the results obtained at 
P > 0.2 GPa (see above and Fig. 2), but the solid was 
unstable in our simulations at P < 0.2 GPa. 

The temperature dependence of the bulk modulus is 
displayed in Fig. 3 for three different pressures. Here 
again open and filled symbols indicate hep and fee he- 
lium, respectively. For each pressure under consideration, 
the plotted data correspond to the temperature region 
where we found the solids to be (meta)stable along the 
PIMC simulations. For each pressure, the bulk modulus 



decreases as the temperature is raised, and this decrease 
is similar for the different pressures. In fact, the obtained 
curves B{T) are parallel within the statistical error of our 
simulations. In the whole accessible temperature range, 
the bulk modulus decreases by 2.1, 2.6, and 3.1 GPa, for 
P = 4, 7, and 10 GPa, respectively. These values amount 
to 13.0, 9.8, and 8.5 % of the corresponding bulk modulus 
at 25 K. 

To connect with data derived from experiment, we note 
that Zha et al. 48] have obtained the bulk modulus of 
solid ^He from Brillouin scattering measurements. From 
these results, they derived the dependence of the isother- 
mal bulk modulus upon material density at 300 K. In 
Fig. 4 we present results of our PIMC simulations for 
^He at 300 K (open squares) along with those derived 
from Brillouin scattering experiments at room tempera- 
ture (solid line). At this temperature, the density region 
between 0.95 and 1.25 g cm~^ corresponds to a pressure 
between 12 and 29 GPa. At p < 1 g cm""^, our results 
coincide within error bars with those derived from Bril- 
louin scattering, and at higher densities the bulk modulus 
derived from PIMC simulations is somewhat larger. For 
example, at density p — I and 1.25 g/cm'^, Zha et al. |48l | 
obtained B — 42.1 and 79.9 GPa respectively, to be com- 
pared with P = 44 and 86 GPa yielded by our PIMC 
simulations of ''He at the same densities and T = 300 
K. This means that at p = 1.25 g cm^"^, the difference 
amounts to about 7%. 



C. Isotopic effect 

An interesting point in this context is the influence 
of the atomic mass on the sohd compressibihty. This 
isotopic effect can be readily obtained from PIMC sim- 
ulations, since the mass is an input parameter in these 
calculations. Similar isotopic effects have been studied 
earlier for the melting curve [20|, |2l| and molar volume 
|23 . |24| of solid helium. For a given material, lighter iso- 
topes form more compressible solids, as a consequence of 
an increase in the molar volume. This is in fact due to a 
combination of the anharmonicity in the interatomic po- 
tential and zero-point vibrations, which are most impor- 
tant at low temperature. Since these quantum vibrations 
are especially large in the case of helium, due to its low 
mass, the isotopic effect on the bulk modulus is expected 
to be appreciable. 

Thus, one expects the bulk modulus of solid '*He to 
be larger than that of ^He. This is in fact the case, as 
derived from our PIMC simulations. In Fig. 5 we display 
the difference AB = B4 — B^ as a. function of pressure 
at T = 25 K. For example, at 1 GPa and 25 K we find 
B = 4.40 and 4.50 ± 0.01 GPa for ^He and "^He, respec- 
tively, and the difference between both isotopes amounts 
to 2.2%. At the same temperature and 7 GPa, we ob- 
tained B = 26.20 and 26.59 GPa, which translates to a 
relative difference of 1.5%. This relative change in bulk 
modulus associated to the isotopic mass is appreciable, 
and slightly larger than that found for the isotopic ef- 
fect on the molar volume of helium crystals. Thus, at 
a pressure of 1 and 7 GPa and T = 25 K, we find a 
volume difference of 1.8 and 0.9%, when comparing '^He 
and ^He [2J]. In general, given a temperature, the rela- 
tive changes in volume and compressibility decrease for 
increasing pressure, since the solid behaves as if it was 
"more classical" [30, 49], and isotopic effects (of quantum 
origin) are consequently reduced. 



D. Thermal expansion 

At this point, it is worthwhile comparing the temper- 
ature dependence of the bulk modulus obtained at fixed 
volume and fixed pressure, from PIMC simulations in the 
NVT and NPT ensembles, respectively. In Fig. 6 we 
have plotted our results for B, as derived in both cases. 
On one side, the simulations in the isothermal-isobaric 
ensemble were carried out at a pressure P = 1 GPa, and 
the results obtained are represented in Fig. 6 as circles. 
On the other side, simulations in the canonical ensemble 
(at various temperatures) were performed for the volume 
obtained in NPT simulations at T = 25 K and P = 1 
GPa. Results of these NVT simulations are displayed as 
open squares. As expected, B obtained in both types of 
simulations agree with each other (within error bars) at 
T < 25 K. At higher temperatures, the NPT simulations 
yield values of the bulk modulus smaller than those de- 
rived from NVT simulations. This is due to the thermal 
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FIG. 5: Isotopic effect on the isothermal bulk modulus as a 
function of pressure. Shown is the difference AB — B4 — B3 
between the bulk modulus of ^He and ''He at 25 K. Filled and 
open symbols represent data for fee and hep helium, respec- 
tively. 
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expansion, that is taken into account in the simulations 
at constant pressure, but is neglected in the constant- 
volume simulations. In general, an increase in crystal 
volume causes a decrease in the bulk modulus. 

In this context, another application of our PIMC sim- 
ulations consists in the calculation of the thermal expan- 
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sion coefficient 
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which can be obtained in the canonical ensemble by using 
the expression ^50i] : 
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This equation in fact relates the thermal expansion with 
the compressibility through the temperature derivative 
of the pressure at constant volume. In connection with 
Eq. ([5]), we note that the pressure obtained in NVT 
simulations (with a fixed lattice parameter a — 3.6955 
A) increases from 1 GPa to 1.3 GPa when temperature 
rises from 25 to 75 K. At 25 K, we have B = 4.46 GPa 
and {dP/dT)v = 8.87 x IQ-* GPa K^^, giving unvt = 
1.99 X lO"** K~^. This result is in agreement with that 
found directly from the volume change obtained in the 
NPT simulations, which yields unpt = 2.03 x 10-^K"^ 
Finally, in Fig. 7 we present results of the thermal 
expansion coefficient a yielded by our simulations in the 
isothermal-isobaric ensemble at two (high) pressures. We 
give values for both •^He and ^He in the hep phase. First, 
we observe an important decrease in a as the pressure is 
raised. For ^He at 300 K, we find a reduction in a by 
a factor of five, when increasing P from 12 to 52 GPa. 
Second, we note that the thermal expansion is smaller 
for "^Hc than for '''He. This is easy to understand, since 



the molar volume of solid '^He is larger than that of ''He 
[23 | , and both volumes have to converge one to the other 
in the high-temperature (classical) limit. 



IV. CONCLUSIONS 

Path-integral Monte Carlo simulations in both NVT 
and NPT ensembles have been shown to be well suited to 
analyse the temperature and pressure dependence of the 
bulk modulus of solid helium. In the isothermal-isobaric 
ensemble, B can be derived from the volume fluctuations 
along a simulation run. In the canonical ensemble, a cal- 
culation of B is more elaborate, but can be carried out 
from thermal averages of various intermediate quantities 
obtained in PIMC simulations. Both ensembles NVT 
and NPT give consistent results for the compressibility 
and thermal expansion of •^He and "^He in the whole re- 
gion of temperatures and pressures considered here. 

At a given pressure, the bulk modulus decreases as 
temperature rises. For pressures between 4 and 10 GPa, 
the change in B has been found to be on the order of 
10%, when temperature increases from the low-T limit 
to the melting temperature of the material. 

Solid ^He is more compressible than '^He. At a given 
T, the difference between bulk modulus of both solids 
increases as pressure rises, but the relative difference be- 
tween them decreases. This isotopic effect on the com- 
pressibility of solid helium is appreciable in the range of 
temperatures and pressures studied here. In fact, at 25 
K and pressures on the order of 1 GPa, it amounts to 
about 2%. 

Apart from a precise calculation of the bulk modulus, 
PIMC simulations in the canonical ensemble can be also 
used to study accurately the thermal expansion of the 
solids under consideration at a given pressure. 
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APPENDIX A: PRESSURE AND BULK 
MODULUS IN THE NVT ENSEMBLE 

1. Pressure 

In the path-integral formalism, the canonical partition 
function for N identical quantum particles obeying Boltz- 
mann statistics (no particle exchange) can be written as 
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where Vij (i = l,...,iV;j = 1,...,Q) are Cartesian co- and from the volume derivative of Zjwt one finds: 
ordinates of the N particles at imaginary time j, (3 — 
{ksT)"^, and $(rij, ...,rjvj) is the potential energy. We 
have defined 



K 



1 / mQ \ 



3NQ/2 



and 



iV! \2'K(ih?) 
mQ 



C = 



2/32^2 



(A2) 



(A3) 



1=1 3=1 

where we have defined 



(A4) 



AT 



H{(3, V, {s,,}) = Cl/2/3 ^ ^(s,,,+i - s. 



1=1 j=i 
+ ^J:HV'^'s,,,...,V'/^s^,) (A5) 

^ 3 = 1 

and the dependence of H on /?, V, and the coordinates 
s.jj has been explicitly indicated. 
The pressure is given by 



1 1 dZjyvT 

^ ]3 Znvt dV ' 



/3V 3V^ ' ^ ' 



(A7) 



where Eh refers to the "harmonic" energy: 



N Q N Q 

Eh = C^5:(r,,+i-r,,)2 = C^^^ ^ 5:(s,,+i-s 

i— 1 J — 1 i— 1 J — 1 



The coordinates are subject to the cyclic condition g^^^j 
Yi,Q+i = Til for all particles i — l,...,N. Changing r^ 
to reduced coordinates s,;,- = V~^^^rii, we find: 



. N Q 
Znvt = KV^^^ / H 11 ^^'^ ^^P \~I^H{I3, V, {s,,})] 

^ A 1 A 1 



(A8) 









ay 



(A9) 



2. Bulk modulus 



The isothermal bulk modulus is 

'dP 



B = -V 



dV 



(AlO) 



and substituting Eq. (jA6p for the pressure, we have: 

V d Zf4YT 



B = VPP^ - 



I3Znvt dV^ 



(All) 



(A6) Taking a second volume derivative of Z^yx, we find: 
I 



1 d Zj\jYT 



pz_ 



NVT 



dV^ 



"«"o-i).^f| -«<?)».) 



ay 



f3V^ 
(Eh^v) 



3^2 \ 3 






(A12) 



1 /a$ 



9y 



-H^v) 



Finally, Eqs. (jAlip and (|A12[) give us the bulk modulus 
_B at volume y and temperature T. 



APPENDIX B: COMPRESSIBILITY IN THE NPT 
ENSEMBLE 

We now have the partition function 



JNPT 



dVe-^P'^Z 



NVT 



(Bl) 



8 



To obtain the compressibility, we calculate the pressure 
derivatives of ItiZj^ipt- We find: 



d In Znpt 
dP 

and 

u In Zj\ipT 



1 uZ^pT 
Znpt dP 



-P{V) 



(B2) 



fdZ_ 



NPT 



Qp2 



72 \ rIP 

NPT \ ^ 



1 u Zj\[p'j' 



Znpt dP^ 



= -p'{VY + p'{V') = -l3'a'y, 



(B3) 



where Oy is the mean square fluctuation of the volume. 
Moreover, from Eq. (jB2p : 



d In ZjsipT 

ap2 



-P 



d{V) 
dP 



(B4) 



and from Eqs. (jB3p and (jB4p . we obtain for the isother- 
mal compressibility k: 



1 (d{V)\ ^Pal 
{V}{dPj^ {V) 



(B5) 



This is the thermodynamic relation obtained in general 
in the isothermal-isobaric ensemble [51]. Thus, it can 
be directly used to obtain the compressibility from the 
volume fluctuations in path-integral simulations, without 
any further calculations. 
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